Critical currents for vortex defect motion in superconducting arrays 
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We study numerically the motion of vortices in two-dimensional arrays of resistively shunted 
Josephson junctions. An extra vortex is created in the ground states by introducing novel boundary 
conditions and made mobile by applying external currents. We then measure critical currents and 
the corresponding pinning energy barriers to vortex motion, which in the unfrustrated case agree well 
with previous theoretical and experimental findings. In the fully frustrated case our results also give 
good agreement with experimental ones, in sharp contrast with the existing theoretical prediction. 
A physical explanation is provided in relation with the vortex motion observed in simulations. 
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In equilibrium the two-dimensional Josephson junction 
array (JJA) is well known to exhibit the Berezinskii- 
Kosterlitz-Thouless transition, 1 which is driven by un- 
binding of vortex-antivortex pairs. Unbound vortices, 
while in motion, dissipate energy and thus make the sys- 
tem resistive. Without an external magnetic field this 
happens above the critical temperature. It is possible, 
however, to induce an unpaired vortex (defect) at low 
temperatures, e.g., by applying weak magnetic fields. 
Such a field-induced vortex defect sits on a plaquette 
where energy takes its local minimum; this is separated 
by the potential barrier set by the underlying lattice. The 
defect may be driven into motion by external currents, as 
it is exerted by the force F — (/i/2e)J x z with J being the 
external current density on the xy plane. When the driv- 
ing force is sufficiently strong, the defect can overcome 
the barrier to move onto the adjacent plaquette and ac- 
cordingly give rise to non-vanishing voltage. Theoretical 
calculation of the barrier height is based on the iteration 
method, static in natureiSi 3 . On the other hand, in exper- 
iment, the barrier height is determined dynamically, by 
measuring critical depinning currents^ 

In the absence of background frustration / = 0, where 
/ measures the number of flux quanta per plaquette, this 
barrier has been estimated as £ ; b(/=0) 0.199 in units 
of the Josephson coupling energy Ejm& The critical cur- 
rent in this case has been obtained from the onset of 
instability of the phase configuration and found to be 
I c (f=Q) w 0.105 in units of the single-junction criti- 
cal current i c , which agrees well with experiment^ For 
high magnetic fields, say, / = 1/2, however, the numer- 
ically estimated value does not agree with the experi- 
mental one: Eb(1/2)/Eb(0) ~ 6 for the former whereas 
E B (l/2)/E B (0) « 1.3 for the latter. This discrepancy 
has remained unexplained since then. 

This work is to bridge the gap between theory and 
experiment: We perform extensive dynamic simulations 
of the resistively-shunted junction (RSJ) model, which 
conveniently describes real dynamics of JJAs. Unlike ex- 
isting studies, which were limited mostly to the defect- 
free case because of the difficulty to create a single defect 



in the presence of external currents, we introduce novel 
boundary conditions for a single extra vortex and study 
its motion in the ground state of the / = 1/2 system as 
well as the / = one. We then measure the critical cur- 
rent I c and the corresponding energy barrier Eb to vortex 
defect motion; this is, to the best of our knowledge, the 
first direct numerical measurement of the barrier height. 
It is confirmed that our results for / = agree well with 
existing numerical and experimental values)^ also re- 
vealed in the fully- frustrated (/ = 1/2) JJA is excellent 
agreement of our results with experimental ones unex- 
plained so far. Via direct observation of actual vortex 
motion, we suggest that the previous overestimation of 
the energy barrier can be attributed to the unrealistic 
configuration of vortices. 

We begin with an L x L square array of RSJs be- 
tween superconducting grains, with shunt resistance R. 
The ith grain, located at position r,-, is described by its 
phase {4>i}- In order to compute the energy barrier, it is 
necessary to create just one extra vortex in the system. 
Both periodic boundary conditions (PBC) and fluctuat- 
ing twist boundary conditions (FTBC) pi used most fre- 
quently, allow / = n/L with an integer n and thus pro- 
duce at the least L extra vortices. This problem can be 
remedied by using the diagonally antiperiodic boundary 
conditions (DAPBC) that introduce 2-7T rotations of the 
phase variables around the whole system and thus pro- 
duce a single extra vortexi However, the straightforward 
application of DAPBC to dynamic simulations fails in 
the presence of external currents since the global current 
conservation cannot be fulfilled^ We thus modify FTBC 
in such a way that the system contains one additional 
vortex. In FTBC, the gauge-invariant phase difference 
is decomposed into two parts: One is periodic across the 
system and thus does not carry voltage whereas the other, 
called the fluctuating twist variable, is directly related to 
the voltage drop across the system. On the former part 
the boundary conditions of the DAPBC type, detailed 
in Ref. 0, are then imposed, relating only the bound- 
ary phase variables on the same columns and rows; this 
makes the application of the FTBC straightforward. 
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Equations of motion for the system (without capacitive 
and inductive couplings) then read, after scaling the time 
in units of h/2ei c R, 
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where the primed summation runs over the nearest neigh- 
bors of grain i, = — r 3 is a unit vector with the 
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lattice constant set to unity, and (f>ij = fa — (, 
together with the voltage-carrying part, describes the 
gauge- invariant phase difference between sites i and j. 
The bond angle A^, given by the line integral of the 
vector potential, takes in the Landau gauge the values 
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with r j denoting the position of site i. The time evolution 
of the twist variable A = (A x , A y ) is governed by£ 
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where denotes the summation over all nearest 

neighboring pairs in the a (= x, y) direction and the dc 
current Id c (measured in units of i c ) is injected along the 
x direction. 

We use the second-order Euler algorithm to integrate 
Eqs. and (J2J with the time step of size At = 0.05. 
As initial conditions, we use the well-known ground- 
state phase configurations obtained with the conven- 
tional PBC: the ferromagnetic configuration fa = for 
/ = and the checker-board configuration faj = ir/4 
for / = 1/2. After appropriate equilibration, the zero- 
temperature vortex configurations in Fig.^are obtained, 
which manifestly show the location of the extra vortex 
and justify the boundary conditions employed. We then 
apply the external dc current Idc and measure the voltage 
drop (in units of i c R) 
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with (• • • ) denoting the time average and the energy 

E - > 'cosifaj - r i:i ■ A). (4) 
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in units of the Josephson coupling strength Ej = M c /2e. 
Typically, data for physical quantities have been averaged 
over 5 x 10 5 time steps after equilibration over 3 x 10 4 
time steps for the JJA with size up to L = 120. 

Figure[21shows the current-voltage (/V) characteristics 
for (a) / = and (b) / = 1/2 in the system of size L = 




FIG. 1: Zero-temperature vortex configurations under the 
boundary conditions in Ref. for / = (left) and / = 1/2 
(right). Filled squares denote the plaquettes occupied by vor- 
tices. 



24, 32, and 64 (from top to bottom), ft is observed that 
nonzero voltage drop develops as Idc is increased beyond 
the critical current I c , the size dependence of which are 
shown in the insets. As L — > oo, the critical currents 
for defect motion unambiguously approach well-defined 
values from aboveS and the two-parameter least-square 
fits yield 

7 c (/=0) = 0.101(1) 
2c(/=l/2) = 0.091(1), 

respectively, where the numbers in the parenthesis denote 
the numerical errors in the last digits. Our obtained value 
7 C (0) is essentially the same as the known valued con- 
firming the validity of our modified FTBC. For / = 1/2, 
however, our result sharply disagrees with the numer- 
ical value I c (l/2) w 0.32 in Ref. and appears much 
smaller than that of the defect-free system^ This indi- 
cates that the critical current to make the point defect 
mobile is much less than the corresponding current for 
the rigid motion of the vortex superlattice. Remarkably, 
our results of the critical current are consistent with the 
experimental result I c w 0.1, below which the dynam- 
ical resistance is negligible for any value of /. 5 It is to 
be noted here that our results are obtained from the di- 
rect numerical calculation of the IV characteristics while 
other studies present mostly indirect estimates from the 
pinning energy barrier with the phase configuration of 
given vorticities assumed (see below) A 

Our simulations also allow us to directly measure the 
pinning energy barrier. To this end, we first plot in Fig.[3| 
the energy in Eq. Q as a function of time in the array of 
size L — 120 at the external current slightly higher than 
the critical value: Idc — I c {f,L) + 0.0001. Observed are 
oscillatory behaviors of E(t) in two very different time 
scales. While the long-period oscillations (only one pe- 
riod of which is shown for simplicity) describe the global 
motion of the defect, the small oscillations shown in insets 
disclose the defect motion in short length scales. The for- 
mer arise from the defect motion across the system, i.e., 
the defect eventually arriving at one end of the system 
and then entering again from the other end. The lower 
energy values at both ends of the system reflect bound- 
ary effects, and the asymmetry in the behavior manifests 
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FIG. 2: IV characteristics at (a) / = and (b) / = 1/2 for 
the system size L — 24, 32, and 64, respectively, from top 
to bottom. Idc and (V') are expressed in units of i c and i c R, 
respectively. Insets show critical currents I c as functions of 
the system size L. 



FIG. 3: Energy E(t) versus time t for (a) / = and (b) / = 
1/2 in the system of size L — 120 at injected currents I^c = 
0.1128 for / = and 0.1001 for / = 1/2. For convenience, 
E(t) has been shifted such that E — is the minimum. Insets 
display expansion plots of E(t) near maxima. 



that the symmetry of the potential energy around the 
middle of the system is broken by the driving currents. 

On the other hand, the small oscillations in the insets 
are directly related with the lattice pinning effects. The 
local minimum (dip) of the energy corresponds to the 
(extra) vortex sitting down at the center of a plaquette 
whereas the peak is produced by the vortex on a link in 
the course of moving on to the adjacent plaquette. In the 
case of / = 0, there appears only one peak in a period. 
Since there are L lattice pinning barriers for a defect to 
move across the system, the longer oscillation in Fig.EJa) 
consists of L small oscillations. For / = 1/2, although 
overall features are similar to those for / = 0, the shape 
of the pinning barrier is markedly different [see the in- 
set of Fig. 13(b)] : One small oscillation consists of double 
peaks separated by a small dip and these peaks are sep- 
arated by bigger dips. The number of bigger dips turns 
out to be L/2, which reveals that the defect motion for 
/ = 1/2 is periodic in space with the period twice the 
lattice constant; this reflects the 2x2 translational sym- 
metry of the vortex superlattice structure in the ground 
state. We define the pinning energy barrier Eb to be the 
energy difference between the bigger dip and the peak, 
neglecting the smaller one. 

To minimize the boundary effects, we measure the 
barrier at the maximum of the long-period oscillation, 



namely when the defect resides near the center of the 
system. Figure 0] presents the measured barrier height 
versus the system size. The extrapolation to the thermo- 
dynamic limit leads toii 

E B (f=0) = 0.192(1) 
E B {f=l/2) = 0.255(1). 

The resulting ratio, Eb(1/2)/Eb(0) ~ 1-33, is again in 
good agreement with the experimental result in Ref. 4. 
Note the dramatic improvement over the existing value, 
E B (l/2)/E B (Q) w 6, obtained via the iteration method. 

We now provide a physical explanation as to the origin 
of the previous discrepancy, in relation with the vortex 
motion directly observed in our simulations. We display 
in Fig. |5] the vortex motion in the / = 1/2 background, 
which reveals the two-step process in the defect motion: 
The defect first pushes away the nearby vortex to make 
the plaquette vacant and subsequently occupies this va- 
cancy. This type of motion gives rise to small oscilla- 
tions of the energy in Fig. [3J In the numerical estimate 
in Ref. |j, however, the vortex defect was assumed to 
sit on the top of the other (background) vortex, which 
results in much higher energy due to strong on-site re- 
pulsion. This type of vortex configuration implies that 
the extra vortex moves in the rigid vortex lattice back- 
ground, which is different from the actual motion of the 
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vortices shown in Fig. [3] We believe that the latter type 
of actual motion occurs in experiment. Namely, vortex 
motion is influenced by background vortices present in 
the frustrated system, thus yielding substantially lower 
critical currents. 



FIG. 4: Pinning barrier height Eb(/) as a function of the 
system size L for (a) / = and (b) / = 1/2. Systems of size 
up to L = 120 are displayed. 



In summary, we have introduced novel boundary con- 
ditions which are convenient to study vortex dynamics 
in JJAs. From dynamic simulations in the presence of 
external currents and magnetic fields, the critical cur- 
rents as well as the pinning energy barriers to vortex 
motion have been directly measured without any a pri- 
ori assumption on the vortex configuration. For / = 0, 
our results agree fully with previous theoretical and ex- 
perimental results. For / = 1/2, our results also give 
excellent agreement with experimental ones, in contrast 
with the previous theoretical prediction. ^From the di- 
rect observation of the vortex motion, it has been found 
that for / = 1/2 the defect moves in the two-step pro- 
cess with two vortices involved, producing double peaks 
in the energy barrier. This suggests that vortices move 
interactively when a large number of vortices are present. 
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FIG. 5: Pattern of the defect motion in the system with / = 
1/2. Currents are applied along the —x direction (i.e., from 
right to left). 
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